close all;
clearvars;
clc;

%% 参数设置
numSubcarriers = 1632;    % 总子载波数（符合NR标准）
numRxAntennas  = 2;       % 接收天线数
numTxPorts     = 2;       % 发送端口数
dmrsType       = 1;       % DM-RS类型（1或2）
symbolIndex    = 1;       % 符号索引

%% 生成DM-RS基序列（Type1示例，间隔2子载波）
X_DMRS = exp(1i*2*pi*rand(numSubcarriers,1)); % 随机基序列（实际可替换为3GPP标准序列）

%% 生成多径信道模型（时域冲激响应）
delaySpread = 10;         % 最大时延采样点数
h_td = (randn(delaySpread, numRxAntennas, numTxPorts) + ... % 时域信道冲激响应
        1i*randn(delaySpread, numRxAntennas, numTxPorts)) / sqrt(2);
    
% 转换到频域
H_true = fft(h_td, numSubcarriers, 1); % [numSubcarriers x numRxAntennas x numTxPorts]

%% 生成接收信号Y_RS（频域）
Y_RS = sum(H_true .* reshape(X_DMRS, [], 1, 1), 3); % 公式：Y = sum(H * X) + noise

% 添加噪声（可选）
SNR_dB = 30; % 信噪比
Y_RS = awgn(Y_RS, SNR_dB, 'measured');

%% 信道估计（全子载波处理）
H_hat = zeros(numSubcarriers, numRxAntennas, numTxPorts);
cdmGroups = floor(numSubcarriers / (dmrsType+1)); % Type1每组占2子载波

for group = 1:cdmGroups
    % Type1子载波对索引（k0, k2）
    k0 = 2*(group-1) + 1;
    k2 = k0 + 2;
    if k2 > numSubcarriers, break; end % 边界保护
    
    % 对每个接收天线处理
    for r = 1:numRxAntennas
        % 计算两个端口的信道估计
        H_hat(k0, r, 1) = (Y_RS(k0,r)/X_DMRS(k0) + Y_RS(k2,r)/X_DMRS(k2)) / 2;
        H_hat(k0, r, 2) = (Y_RS(k0,r)/X_DMRS(k0) - Y_RS(k2,r)/X_DMRS(k2)) / 2;
        
        % 保持组内子载波信道估计相同
        H_hat(k2, r, 1) = H_hat(k0, r, 1);
        H_hat(k2, r, 2) = H_hat(k0, r, 2);
    end
end

%% 时域转换（取端口1天线1的估计结果示例）
h_est_td = ifft(H_hat(:,1,1)); % 时域冲激响应估计



%% 绘图
figure('Position', [100 100 1200 600])

H_true_used = H_true(1:2:end,:,:);
H_hat_used = H_hat(1:2:end,:,:);

% 频域幅度响应
subplot(2,2,1)
plot(abs(H_true_used(:,1,1))), hold on
plot(abs(H_hat_used(:,1,1))), hold off
title('频域幅度响应（端口1天线1）')
xlabel('子载波索引'), ylabel('|H(f)|')
legend('真实信道','估计信道')

% 频域相位响应
subplot(2,2,2)
plot(angle(H_true_used(:,1,1))), hold on
plot(angle(H_hat_used(:,1,1))), hold off
title('频域相位响应（端口1天线1）')
xlabel('子载波索引'), ylabel('∠H(f)')

H_td_used = h_td(1:2:end,:,:);
H_est_td_used = h_est_td(1:2:end,:,:);

% 时域冲激响应幅度
subplot(2,2,3)
stem(abs(H_td_used(:,1,1)), 'filled'), hold on
stem(abs(H_est_td_used(1:delaySpread)), '--'), hold off
title('时域冲激响应幅度')
xlabel('采样点'), ylabel('|h(t)|')
legend('真实信道','估计信道')

% 时域误差
subplot(2,2,4)
plot(abs(H_td_used(:,1,1) - H_est_td_used(1:2:delaySpread)))
title('时域估计误差')
xlabel('采样点'), ylabel('绝对误差')

%% 性能评估
mse_freq = mean(abs(H_true_used(:) - H_hat_used(:)).^2);
mse_time = mean(abs(H_td_used(:,1,1) - H_est_td_used(1:2:delaySpread)).^2);
disp(['频域MSE: ', num2str(mse_freq)]);
disp(['时域MSE: ', num2str(mse_time)]);